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I. INTRODUCTION 

Experiments of Bose-Einstein condensation (BEC) 
with ^Li atoms is still a challenging research area as the 
s-wave scattering length (os) is negative, which indicates 
an attractive atom-atom interaction [1] . A homogeneous 
condensate with negative scattering length is impossible, 
as it leads to ever increasing density [2]. As Ug is negative, 
the attractive interaction energy gradually increases with 
increase in the number of atoms in a small volume and 
the condensate approaches collapse. However, the situ- 
ation changes drastically in the presence of an external 
confinement. A spatially confined BEC may be stable for 
a small, finite number of atoms [3] . In the presence of a 
confining potential, the destabilizing effect of the effective 
negative interaction is balanced by the kinetic pressure 
of the gas and a metastable condensate can form [1] . For 
■^Li, as = -(27.3 ± 0.8) Bohr and for T=0, a metastable 
condensate exists when the number of atoms is less than 
the critical number, which is roughly 1300 [4], whereas 
theories predict that BEC can occur in a trap with no 
more than about 1400 atoms [3,5-7]. 

During the last decade a number of articles have been 
published where the properties and stability of the at- 
tractive condensate have been discussed in details [3,5-7]. 
The two- and three-body decay processes near the col- 
lapse have also been studied in a number of papers [7,8]. 
However, all the earlier calculations use the mean-field 
approximation and the ground state energy is calculated 
by the Gross-Pitaevskii (GP) energy functional [^. GP 
theory is based on the pseudopotcntial form of the atom- 
atom interaction i.e. a (5- function potential, where the 
strength of interatomic interaction is absolutely deter- 
mined by a single parameter a^. But due to the presence 
of the pathological singularity at \x\ = 0, the Hamilto- 
nian becomes unbound from below and it has been em- 
phasized previously that a J-function is not suitable as an 
exact potential in 3D attractive systems [3, [io| ■ Again, 
as the attractive BEC becomes highly correlated near the 
critical point, uncorrclated GP equation can not take care 



of the effect of interatomic correlation. 

Thus the motivation of the present work is to study 
the same system using an ab initio many-body method, 
that is the potential harmonic expansion method |lll.ll2l|. 
with the incorporation of a realistic potential, like the 
van der Waals potential. The presence of a hard sphere 
below some cutoff radius and the —Ce/\x\^ tail prop- 
erly take care of the effects of both short-range and 
long-range correlation and accurately represent the real- 
istic features. We find the commonly known low-density 
metastable branch and discuss its stability, structure and 
decay rates due to inelastic two- and three-body colli- 
sions. Besides the known solution, we also find another 
branch of solution at a higher density. Due to the use 
of the realistic van der Waals interaction in the many- 
body calculation, we have a deep attractive well on the 
left side of the metastable region in the effective poten- 
tial. In the standard GP theory, the metastable region 
just vanishes at the critical point and the whole conden- 
sate collapses into the singular well and the fate of the 
condensate is not predicted further. However, in our cal- 
culation the metastable condensate leaks through the in- 
termediate barrier and settles down in the extremely nar- 
row (width ^ O.Obfim) well. The atoms become highly 
correlated and due to high two-body and three-body col- 
lisions, atoms form a cluster. This corresponds to the 
high-density branch. In the present communication, we 
investigate the transition between the two branches of 
solutions as a function of the number of atoms and their 
dependence on the trap size. We also calculate decay 
rates of the condensate due to two- and three-body col- 
lisions. In other articles [1^, the existence of a similar 
type of new stable branch has been reported, which is 
intermediate in density between the dilute metastable 
state and the collapsed state. This is described as the 
effect of use of non-local interaction in the GP theory. It 
has been pointed out that for ^Li system, the scattering 
cross-section has a momentum dependence and the effec- 
tive interaction is non-local, changing from attractive to 
repulsive at a characteristic range r^. They calculated 



the properties of the attractive condensate by the vari- 
ational technique. Using a Gaussian trial wavefunction 
they minimized the quantum action. 

In the present communication, we also compare 
our results with those obtained using the non-local 
interaction. A road map of the present study is given 
below. Sec. II contains the many-body approach used 
in this work, which is based on the potential harmonic 
expansion method. Sec. Ill presents the calculation with 
a non-local interaction. Sec. IV contains a comparison 
of the structure and stability of the bosonic system of 
both the low- and high-density branches obtained from 
two different theoretical calculations and the discussion 
of the many-body effect. Sec.V provides a summary of 
our results and their conclusions. 



II. MANY-BODY CALCULATION: (A) 
POTENTIAL HARMONICS BASIS 

We consider a system of A=(N+1) identical spinlcss 
bosons (each of mass m) confined in an external trap- 
ping potential Utrap{xj) and interacting via a two-body 
potential Y{xi — Xj), Xi being the position vector of the 
i-th particle. The corresponding Schrodinger equation is 
written as 



2m 



5I^' + I]^*™p(^0 



i=l 



4=1 



^ V{x, ~ X,) ~ E\ij{xi, ...,xa) = Q. (1) 



We define N Jacobi vectors as 



2i 



Xi+l 



-S^) 



(z = l,...,iV). (2) 



The center of mass is ^ = Si=i ^i/-^- Since the labelling 
of particles is arbitrary, we choose the relative separation 
of {i,j) interacting pair {xij ~ Xi ^ Xj) as Cat. We next 
define the hyperradius r of the set of N Jacobi vectors 
through [m 



N 



^^ = Ecf = 7E4-2E-? , 



(3) 



i,j>i 



where r^ = Xi — R is the position vector of the i-th particle 
from the center of mass of the system. In this way the 
relative motion of the bosons is described by 



N 



E^c. +^t™pW + 



V^nt{Cu-,CN)-En ^(Ci,...,Cw)=0, (4) 



where En is the energy of the relative motion. In the 
present work Utrap{xi) = ^muj'^\xi\'^ is a spherically sym- 
metric harmonic oscillator potential and consequently 
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Utrap{xi) = -muj r + -mAuj R 



(5) 



The first term on the right side is the effective trap po- 
tential Vtrapi'f) for the relative motion and the second 
term is the trap potential of the center of mass motion. 
The equation for the center of mass motion separates 
completely and is simply the equation for a three dimen- 
sional harmonic oscillator. The total energy of the sys- 
tem is thus {Er + ^hw). In Eq. (4), Vint is the sum of 
all pair-wise interactions, expressed in terms of the Ja- 
cobi vectors. In our approach we consider that when [ij) 
pair interacts, the rest of the bosons are inert spectators. 
So we define a hyperradi us / o,;, for the remaining (N-1) 
noninteracting bosons as [13 | 
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(6) 



so that r^ = x 



The hyperangle (f> is introduced 



through Xij = r cos (j) and pij = r sin (f>. Besides r, 0, i9, (p 
(where d and Lp are the polar angles of x^), there are 
(3iV — 4) remaining variables. These are constituted 
by 2{N — 1) polar angles associated with (^i, ..., ^Af-i) 
and [N — 2) angles defining their relative lengths and 

collectively denoted by ^)v_i' called hyperangles in the 
3(iV — l)-dimensional space. The corresponding form of 
Laplace operator can be found in [ll[. In the potential 
harmonics expansion method (PHEM), -0 is decomposed 
into Faddeev components, (pij for the {ij) interacting pair 
as, 

A 

V'= E '^•y- 

Then Eq. (4) can be written as 



(7) 



(T + Vtrap - Eji)4)ij == -V{xij) ^ 0A;;, 



(8) 
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where T 



"■^X]i=i^^- Note that the assumptions 



that correlations higher than two-body ones in ^/j are 
negligible and that the angular and hyperangular mo- 
menta of the system are contributed by the interacting 
pair only [ig , make the Faddeev component 0ij indepen- 
dent of the coordinates of all the particles, other than the 
interacting pair |ll| . Thus (f)ij = 4>ij {xij , r) . With this 
assumption we expand (f>ij in the subset of hyperspherical 
harmonics (HH) necessary for the expansion of V{xij). 
Thus, 



4>ij{x.i.j,r) = r 



-(3iV-l) v^^/' 



K 



(9) 



This subset of HH is called the potential harmonics (PH) 

basis and is denoted by {'P2k+i{^^n)}'' ^jv denotes the 
set of hyperangles in 3N dimensional space with the inter- 
acting pair uniquely identified as above. Note also that 

'^2K+ii^%) i^ independent of (Ci, ■■■, Ca^-i) and thus no 
contribution to the orbital and hyper angular momenta of 
the condensate conies from the remaining (TV — 1) non- 
interacting bosons. Thus the total orbital angular mo- 
menta of the condensate and its projection are simply 
those of the interacting pair, viz. I and m respectively. 
The complete analytic expression for 'P2'ff-i-/(^w) '^^n be 
found in [IJ] . Substitution of Eq. (9) into Eq. (8) and 
subsequent projection on a particular PH gives [la] 



h"^ <P h^ Ck{Ck + 1) 
m dr^ m r^ 



+ Vtrap{r)- En^ u^(r) 



+ Y.fK'iVKK'{rWj,,{r)^0, 



(10) 



K' 



where Vkk' {i^) is the potential matrix and is given by, 



Vkk' (r) 



■ylrn* 
2K+1 



{n'^)v{x,,)r^^,^,in'i)dn'^. (ii) 



The quantities Ck and /^j are given by 



Jki 



■iN-Z 



= 2K + 1+ 2 
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2A'+/(^V)k 2K+z(^^7v) >' 



(12) 



the latter being the overlap of the PH for the (zj)- 
partition (corresponding to only the (ij)-pair interacting) 
with the sum of PH for all partitions. The quantum num- 
ber K in Eq. (10) is the hyperspherical (grand orbital) 
quantum number in the 3 A^ dimensional space |lll . Il2| , 
analogous with the orbital angular momentum in three 
dimensional space. The ordinary orbital angular momen- 
tum of the system and its projection are I and to, as indi- 
cated earlier. All other intermediate angular momentum 
quantum numbers of the system take zero eigen values, 
due to the choice of the PH basis. Thus the choice of the 
PH basis reduces the complications of the many-body 
system immensely. Eq. (10) can be put in a symmetric 
form as 



TO dr"^ 



;{C{C + 1) + AK{K + a + P+l)} 



+ Vtrap{r)-ER\^UKl{r) 



K' 



where C = Z + (3A-6)/2, UKi{r) = fKiu'j,{r) 
and/3 = l + h- 



(13) 
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III. MANY-BODY CALCULATION: (B) 

INCORPORATION OF TWO-BODY 

CORRELATION FUNCTION 

In the experimentally achievable BEC, the average in- 
terparticle separation is much larger than the range of 



two-body interaction. This is indeed required to prevent 
atomic three-body collisions and formation of molecules. 
Moreover, the energy of the interacting pair is extremely 
small. Thus the two-body interaction is generally repre- 
sented by the s-wave scattering length (a,,). In the GP 
picture, the effective interaction is given in terms of as 
alone, making the results independent of the shape of 
the interatomic potential. Hence, a positive value of Qs 
gives a repulsive condensate and a negative Os gives an 
attractive condensate. However, a realistic interatomic 
interaction, like the van der Waals potential, is always 



associated with an attractive — 1/ 



tail at larger sep- 



arations and a strong repulsion at short separations. De- 
pending on the nature of these two parts, Os can be either 
positive or negative |14{ . For a given two-body potential 
V{xij) having a finite range, a^ can be obtained by solv- 
ing the zero-energy two-body Schrodinger equation for 
the wave function ri{xij) 



ft2 1 d 

TO xjj dXij 



^2dv{xij) \+v{x,Mx,,) = 0. (14) 



dxi 



The asymptotic part of r]{xij) has the form (ciXy -I- C2) 
and the corresponding s-wave scattering length is given 
by a, = -fa. 

Now the rate of convergence of the PH expansion in 
Eq. (9) is very slow. This can be understood from the 
fact that the leading term of this basis, viz. the term with 
ii" = 0, is a non-zero constant, while 0ij must be vanish- 
ingly small for small values of Xij due to the strong short- 
range repulsion of V{xij). Consequently a large number 
of terms is needed to represent 0y faithfully for small val- 
ues of Xij . The rate of convergence is improved dramati- 
cally by introducing a short-range correlation function in 
the expansion basis [l^, that represents the true nature 
of (f>ij{xij) as Xij — > 0. In the experimental BEC, the en- 
ergy of the interacting pair (which is ~ hjj) is negligible 
compared with the depth of interaction potential V{xij) 
(which is of the order of typical atomic energy scale). 
Thus Tjixij) is a good representation of the short-range 
behavior of (j)ij [xij ) . Hence we use 'q{xij ) as a short-range 
correlation function in the PH expansion, to enhance its 
rate of convergence [13] : 



4>ijixij,r) =r~ 



E^: 



2K+ii^%WKirHxv] 



(15) 



K 



As ii{xij) correctly reproduces the short separation be- 
havior of the interacting-pair Faddcev component, con- 
vergence rate of the PH expansion for the residual part 
of ipij is quite fast. Validity of this statement has been 
tested in numerical calculations. In an actual calculation, 
we solve Eq. (14) for the zero energy two-body wave func- 
tion ri{xij) in the chosen two-body potential V{xij), after 
adjusting the hard core radius (jc), such that Oj has the 
desired value [l^. This ri{xij) is then used in Eq. (15). 



Then the potential matrix becomes 



1 + 2 



P^^{z)v(r^m^jwi{z)}dz. (16) 

where P^ (z) is the Jacobi polynomial, whose norm and 
weight function are h'^ and wiiz) respectively [3. We 
restrict the K-sum in Eq. (15) to an upper limit Kmax 
from the requirement of convergence. Finally we solve 
the set of coupled differential equations (CDE), Eq. (13), 
by hyperspherical adiabatic approximation [20|. The to- 
tal energy of the system is obtained by adding energy of 
the center of mass motion {3haj/2) to Ej^. The wave func- 
tion of the condensate in the physical space can be ob- 
tained from UKi{r), using Eqs. (7) and (9) and the rala- 
tions between hyperspherical and physical coordinates. 



IV. GP THEORY WITH NON-LOCAL 
INTERACTION 

In the standard treatment of alkali-metal atoms, the 
interatomic interaction is chosen as a local form, i.e., 
momentum independent. The effective zero-range pseu- 
dopotcntial has the form V{x) = g6^{x), where g = 
^'^^ °° . However, this is not correct when the scattering 
cross-section has a significant momentum dependence. 
This implies that the effective interaction is non-local, 
changing from attractive to repulsive at a characteristic 
range r^. The system of ^Li atoms exhibit such momen- 
tum dependence and the effective ^Li-^Li interaction can 
be written as [11] 



V{k) 



Airh^ 



[ar + {as - ar)f{kre)] 



(17) 



where Cg corresponds to the attractive potential and has 
the value a^ = —27.3 Bohr, whereas the repulsive contri- 
bution is modeled by a local positive term with a,- = 6.6 
Bohr and Ve = 10^ Bohr. The shape function f{x) may 
be chosen as a Lorenzian f{x) = {l + x'^)~'^. In real space 
the effective inter-atomic potential then reads 



V{x) = -6^{x) H ^ 



a^) e' 



\x\/rc 



(18) 



This potential appears in the nonlocal GP energy func- 
tional 
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rix) 



2m 



Utrapi\x\) 



{x) d^x 



|(/)(f)|V(f - x')|0(f' )r d^x d^x' , (19) 



where 0(a;) is the spherically-symmetric GP wave func- 
tion. In earlier calculations |13j the variational ansatz of 



the wave function is chosen as 

1 



<P{\x\)=N'^^- 



-3/4^3/2 



exp 



2a2a2 



(20) 



where a is the variational parameter, aao being the width 
of the condensate with an = \/-^ the characteristic 
length of the harmonic trap. In this way the energy E of 
the system becomes [l^ 



Eia) 



Nhjj 



Scr^ + 2N 



lr_ 
3 



t-25(xct) 



(21) 
with "fr = {2/TTy/^ar/ao, n = {2/TT)^^^ao{as - ar)/rl, 
X = 2"^/2ao/re, T2 = al{as - ar)/rl and g{x) = 
exp(a:^)[l — erf(x)], erf(a;) is the error function. The ex- 
trema of E[a) are obtained from the following expression. 



N ={\-a'^) 



-IrO 



no- 2xt2(t3 _ 2x^T20-'^ff(xg) 
3 30F 



(22) 

which actually gives the number of bosons as a function 
of the scaled size a of the cloud. Depending on the choice 
of parameters and the number of atoms, this equation ei- 
ther has one root or three roots. When three solutions are 
present, one corresponds to the low-density metastable 
solution (which we expect for local interaction), one cor- 
responds to the high density stable solution (which is the 
effect of non-locality) and the middle one corresponds to 
an unstable state {i.e. local maximum). The variational 
results for different trap sizes are discussed in the result 
section. 



V. RESULTS 

We consider the experiment on ^Li condensate at Rice 
University pj, for which as = —27.3 Bohr. We take oq as 
the unit of length and energy is expressed in units of the 
oscillator energy quanta huj. The trap size correspond- 
ing to the Rice University experiment is ag = 3.0fim. 
For the many-body calculation, we choose the van der 
Waals potential, whose short range repulsion is charac- 
terized by a hard core of radius r^ and the long range 
part is described by an attractive tail —Ce/\xij\^, where 
Cg is the strength parameter which is known from exper- 
iments. Choice of a realistic interatomic interaction like 
the van der Waals potential is very important, as it has 
been shown that the shape-independent approximation 
is not valid in typical laboratory BECs [21|. In the local 
mean-field description the two-body interaction is usually 
represented by the s-wave scattering length a^. To mimic 
the ^Li trap of Rice University, our chosen parameters are 
Cg = 1.71487 X 10-12 o.u. and r^ = 5.3378 x 10"^ o.u., 
which reproduce the experimental value of a^. With these 
sets of parameters we solve the set of coupled differen- 
tial equations by the hyperspherical adiabatic approx- 
imation (HAA) [201. We assume that the hyperradial 



motion is slow in comparison with the hypcrangular mo- 
tion. Hence, for a fixed value of r, the equation for the 
hyperangular motion can be solved adiabatically. The en- 
ergy eigen value of this equation is a parametric function 
of r and provides an effective potential for the hypcrra- 
dial motion. In the HAA prescription, the lowest lying 
such potential is used for the ground state of the sys- 
tem. Although Eq. (13) involves the hyperradius r only, 
the hyperangular motion appears through the coupling 
matrix Vku' ('") • Solution of the hyperangular motion for 
a fixed value of r is equivalent to diagonalizing the hy- 
perangular Hamiltonian in the potential harmonics basis, 
i.e. diagonalizing the potential matrix Vkk'(t) together 
with the hypercentrifugal term of Eq. (13). The lowest 
eigenvalue of this matrix, wq ('') i is the effective potential 
in the hyperradial space, in which collective motion of the 
condensate takes place. Thus the energy and wave func- 
tion of the condensate is obtained by solving the equation 
for hyperradial motion 






K„ 
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' dr 



En 



Co(r) - 



(23) 

where the third term is the overbinding correction 
term [20| and xko is the K-th component of the column 
vector corresponding to the eigenvalue iOo{r). Energy 
and wave function of the system are obtained by solv- 
ing Eq. (23) numerically, subject to appropriate bound- 
ary conditions. The partial wave Uxiir) appearing in 
Eq. (13) is given in HAA as Umir) « Co{r)XKo{r) M- 
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FIG. 1: (Color online) Plot of effective potential cJo('") in o.u. 
against r in o.u. for 1000 ^Li atoms with ao = 3.0/xm. Panel 
(a) shows the narrow attractive well, while panel (b) displays 
the metastable region in detail. Note the different horizontal 
and vertical scales in the two panels. 

In Fig. 1 we plot the effective potential, uJo{r) for 
N — 1000 atoms of ^Li, the condensate being metastable 
as N<Ncr. For N<Ncr, the intermediate metastable re- 
gion (MSR) of ujo{r) is bounded by the high wall of the 
external trap on the right side and on the left by an in- 
termediate barrier (IB) [panel (b)]. On the left of the 
IB, a very deep narrow attractive well (NAW) appears 
which is again followed by a steeply increasing repulsive 
wall, as r decreases [panel (a)] . The existence of a strong 
repulsive core near the origin (r — > 0) is in sharp contrast 
to the GP mean-field picture with local interaction and 
is the immediate reflection of the use of realistic van dcr 



Waals interaction in quantum many-body calculation. In 
the GP picture, negative Og presents a singular attraction 
(with no repulsive core) near the origin. The combination 
of repulsive core and NAW in the present approach, pre- 
vents the Hamiltonian from being unbound from below. 
As the number of atoms increases towards the critical 
number, the height of the intermediate barrier decreases, 
the NAW becomes more negative and narrow. At N = 
Ncr, the MSR along with the intermediate barrier dis- 
appear, with the local maximum of the IB and the local 
minimum of the MSR merging to form a point of inflex- 
ion [l3|. The value of N where this occurs is Ncr- As N 
becomes larger than Ncr, ah the atoms get trapped in the 
NAW, giving rise to a cluster state resulting from the in- 
creased interatomic correlations and two- and three-body 
collisions. The release of binding energy is observed as 
a burst of energy in experiments. This is the collapse 
scenario in the present approach. 

Thus the qualitative features of the pre-collapse sce- 
nario are in fair agreement with the GP local mean- 
field picture away from the critical number. When N 
= Ncr, the energy per particle becomes — oo in the lo- 
cal GP method, and the Hamiltonian becomes unbound 
from below. However, in our many-body calculation the 
Hamiltonian has a lower bound due to the presence of a 
hypercentrifugal barrier as r — >■ along with the use of a 
realistic interatomic interaction with a repulsive core of 
finite size. In the post-collapse state the atoms accumu- 
late in the narrow attractive well and form clusters. 




Q^/1^1o I" I -1000 

0.4 0.8 1.2 1.6 2 20 40 60 80 100 



500 1000 1500 2000 
N 
(a) 




500 1000 1500 2000 
N 
(b) 



FIG. 2: (Color online) Size of the condensate as a function 
of the number of ^Li atoms using PHEM with van der Waals 
interaction and GP with local interaction (having same as) 
[panel (a)], and GP with non-local interaction [panel (b)]. 



We calculate the size of the condensate as a function 
of the number of ^Li atoms, for the choice of the param- 
eters mentioned earlier. We define the average size of 
the condensate (vav) as the root mean square distance 
of individual atoms from the center of mass and is given 
by El 
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< 7-2 >l/2 
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(24) 



where R is the center of mass coordinate. We present 
average size of the condensate calculated by PHEM in 
panel (a) of Fig. 2. The upper branch corresponds to the 



metastable condensate. With increase in particle num- 
ber, the total attractive interaction increases as — ^ ^~ , 
the system contracts and rav decreases as expected. For 
comparison we present in panel (a) the results of the GP 
equation with contact interaction, i.e. with local interac- 
tion. The lower branch, which corresponds to the high- 
density stable branch in the deep attractive well, starts 
from N=20. In the GP mean-field approach with local 
interaction, there is no other stable branch after collapse 
as the whole condensate falls into the singular well. In 
our calculation we find that the size of the condensate in 
the deep attractive well is of the order of 0.003 fim, which 
is the order of the size of the atomic cluster. This clearly 
indicates that the metastable condensate forms clusters 
due to high two-body and three-body collisions within 
an extremely narrow well of width ^ 0.05 fim. However, 
the transition from upper branch to lower branch is dis- 
continuous. In panel (b) of Fig. 2, we present results for 
the same trap but with non-local interaction as described 
in Sec. III. The quantity aao in Eq. (20) is basically the 
average size of the condensate and is denoted by rav in 
panel (b) of Fig. 2. The variational parameter a is cal- 
culated as a function of the number of bosons using the 
algebraic equation [Eq. (22)]. The effective potential en- 
ergy E{a) [Eq. (21)] has the same qualitative features 
as the many-body effective potential u>o{r). Unlike the 
GP effective potential with local interaction, (having one 
local minimum), the use of non-local interaction in the 
calculation of effective potential offers an additional ab- 
solute minimum together with the local minimum. This 
absolute minimum corresponds to the high-density stable 
solution whereas the local minimum corresponds to the 
usual metastable solution. The transition between the 
low-density and high-density branches is again discon- 
tinuous as observed in many-body calculation. However, 
unlike panel (a) of Fig. 2, where the high-density sta- 
ble branch appears from quite small values of N = 20, in 
panel (b), the lower branch starts from N ~ 234. This dis- 
agreement is due to the fact that the origin of non- locality 
and the occurrence of absolute minimum is different in 
the two approaches. In the many-body calculation the 
effect of non-locality is inherent due to the presence of 
hypercentrifugal repulsion (arising from the many-body 
treatment) and the use of van der Waals potential with a 
finite range hard core part whereas in the other approach 
the non-locality is created by assuming that the attrac- 
tive potential has a finite range Vg together with a local 
positive term defined by high-energy scattering length 
Ur > 0. This is reflected by the presence of a non-local 
minimum in the effective interaction. 

Next, to study the nature of discontinuity between the 
two branches of solution, we repeat our calculation for 
other trap sizes. We plot the PHEM results in panels 
(a) - (c ) of Fig. 3 with ao == O.l^im {u = 907.88kHz), 
ao = 0.05/iTO (w ~ 3.63MHz) and ap = O.Ol/im (w = 
90.79MHz). We observe that by reducing the trap size 
the discontinuity is strongly reduced and at ao = O.Olfim, 
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FIG. 3: (Color online) Size of the condensate as a function of 
the number of ^Li atoms using PHEM [panels (a) - (c)] and 
GP with non-local interaction [panels (d) - (f )] . Note that in 
panels (a) - (c), the horizontal axis is upto A'^ = 500, to show 
the small N part prominently. Calculated values of rav for N 
above 500 remain practically constant for PHEM. 



the discontinuity disappears. There is a smooth evolu- 
tion from a dilute cloud to a less dilute cloud. We also 
observe that the size is almost independent of the trap 
size for N > 200, as expected for clusters. For compari- 
son, in panels (d) - (f) of Fig. 3, we plot the variational 
results with non-local interaction for various trap sizes 
ao = 0.6/iTO, ao = 0.3/im and ao = 0.1/j,to. 
We get almost same qualitative pictures as shown in 
panel (a) - (c) of Fig. 3; reducing the trap size, the dis- 
continuity is reduced and at ao = O.I/ztti, the metastable 
branch disappears. However, the quantitative disagree- 
ment between many-body results and variational results 
with non-local interaction is again attributed to the dif- 
ferent origin of non-locality in the effective potential. 
In the variational calculation with non-local interaction, 
the discontinuity disappears at ao — O.lfim whereas the 
corresponding trap size in the many-body calculation is 
ao = O.Olfim. However, in both cases we observe that 
the size of the denser cloud is independent of large N, 
which implies that the atoms are being self-trapped in 
the non-local minimum. 

For the case of a spherical trap with trap size ao = 
3.0/X771, we fail to find the metastable numerical solutions 
beyond N ~ 1430 atoms. The instability occurs as the 
number of atoms increases, the net attractive interaction 



between atoms becomes dominant and kinetic energy can 
no longer stabilize the wave function. In the mean-field 
approach it has already been pointed out that the tran- 
sition to an unstable state may occur due to quantum 
tunnelling. This tunnelling is similar to the ordinary 
quantum tunnelling and its rate can be calculated by 
the semiclassical formula. As this rate is quite small, ne- 
glecting tunnelling we consider only the loss rates due 
to two-body dipolar collisions and three-body recombi- 
nation collisions. 



The upper branch accounts for the loss rate in the deep 
narrow attractive well. As all the atoms are now con- 
fined within a very narrow and deep well, they suffer 
vigorous two- and three-body collisions and the calcu- 
lated loss rate is much higher than the same calculated in 
the mctastable region. The transition between these two 
branches is again discontinuous. In panel (b) of Fig. 4 we 
plot r^ calculated using variational trial wave function 
of Eq. (20) using non-local interaction. Substitution of 
trial wave function Eq. (20) in Eq. (25), leads to 
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FIG. 4: (Color online) Loss rate due to two- and three-body 
collisions as a function of the number of ^Li atoms using 
PHEM [panel (a)] and GP with non-local interaction [panel 
(b)]. 

The total loss rate is given by 

Tl =K f dPx\(l){x)\'^ + L f d^x\<f>{x)\^ (25) 



where two-body dipolar loss rate coefficient K = 1.2 x 
10^^/im^sec~^ and the three-body recombination loss 

10^^/im^sec~ 



rate coefficient L ~ 2.6 x 10^ /im°sec~ [2^. 0(a;) is 
the condensate wave function. This is simply given by 
Eq. (20) for the non-local GP approach. For the many- 
body approach, the wave function ?/; is initially obtained 
as a function of the Jacobi coordinates. It is then trans- 
formed into a function of position vectors {x\,X2, ..., xa} 
of A particles. Finally, the one-body density is obtained 
by integrating ji/jj^ over all Xi, except one. This one- 
body density is used for |(/>(x)|^ in Eq. (25). It is already 
known that in case of a negative scattering length, the 
rapid increase in the density of atoms and large increase 
in the net attractive interaction for a large number of 
atoms, result in a very large loss rate near the critical 
point. Thus for the unique metastable branch in GP 
mean-field picture with local interaction, the loss rate in- 
creases very fast with the increase in the number of atoms 
[7]. However, due to the presence of the deep attractive 
well in the many-body picture, we calculate the loss rate 
for both the branches and observe the discontinuity. By 
using condensate wave function as described above, we 
calculate Fl by using Eq. (25) and the results are plot- 
ted in panel (a) of Fig. 4. The lower branch corresponds 
to the metastable region. As the metastable region is 
flatter compared to the NAW, although the atoms are 
correlated, the probability of two or more atoms to come 
very close is relatively small. Thus the loss rate is ini- 
tially small and increases with increase in atom number. 
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In this formula, N is related to a through Eq. (22). 
The existence of low-density metastable branch and high- 
density stable branch which are separated by an interme- 
diate unstable branch is also reflected in Fig. 4. It is in 
sharp contrast to the existence of a unique metastable 
branch in the mean-fleld approach with local potential. 
For the sake of completeness, in Fig. 5 we plot the loss 
rate for traps of different sizes. 
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FIG. 5: (Color online) Loss rate due to two- and three-body 
collisions as a function of the number of ^Li atoms using 
PHEM [panels (a) - (c)] and GP with non-local interaction 
[panels (d) - (f)]. Note that in panels (a) - (c), the horizontal 
axis is upto A'^ = 500, to show the small A*' part prominently. 
Calculated values of Tl for N above 500 remain practically 
constant for PHEM. 



VI. CONCLUSIONS 

In summary, we have analyzed the stabihty of the at- 
tractive Bose-Einstein condensate of ^Li atoms in a har- 
monic trap using an ab initio quantum many-body cal- 
culation incorporating realistic van der Waals interaction 
and all possible two-body correlations. Unlike the single 
mctastablc condensate arising from the GP mean-field 
theory with local interaction, we observe that the atomic 
cloud may exist in different states. In addition to the 
commonly known low density metastable BEG state, we 
find a high density stable solution, which corresponds to 
the cluster state of the lithium atoms. Between these two 
minima, there is a maximum of the effective potential, 
which does not correspond to any stable physical state 
of the Bose gas. The local GP approach produces only 
the metastable solution, together with an intermediate 
barrier (local maximum), but the effective potential has 
no lower bound on the left of the intermediate barrier. 
Indeed there is a pathological singularity at the origin, 
which only produces a pathological collapse (not leading 
to any physical solution). The effect of a non-local mo- 
mentum dependent interaction (replacing the local con- 
tact interaction) in the GP approach was investigated 
in earlier works |I3j using a variational ansatz. The re- 
sults of this non local GP approach agree qualitatively 
with our present observations. The lack of quantitative 
agreement is attributed to the difference in the origins of 
non-locality in the two approaches. We have investigated 
the average size of the bosonic system, as also the loss 
rates due to two- and three-body collisions for both the 
metastable BEG state and the collapsed cluster state. In 
general these quantities lie on separate branches for the 
two states. For larger trap lengths, the two branches are 
unconnected. But there is a continuous transition for 
very tight traps. The local GP approach produces only 
one branch, corresponding to the metastable state. The 
local GP uses an attractive contact interaction for the ^Li 
condensate. This causes a pathological singularity at the 
origin, preventing any stable collapsed state. The non 
local GP approach uses a non-local potential, which is 
attractive at larger separations, but has a repulsive con- 
tact interaction as r — >■ 0. Hence there is no pathological 
singularity and the effective potential becomes strongly 
repulsive as r — > 0. This gives rise to the stable high den- 



sity branch. In the PHEM approach, the realistic van der 
Waals (vdW) potential is used. This has a very strong 
repulsion at short separations and an attractive I/|xij|^ 
long tail. The strongly repulsive core of the vdW poten- 
tial, together with the hypercentrigugal repulsion, arising 
from the PHEM approach, produces a repulsive core fol- 
lowed by a deep attractive well as one moves away from 
the origin in the effective potential, u)Q{r). Thus both 
the PHEM and the non local GP approaches produce 
the same qualitative results. However, quantitative dif- 
ferences exist due to different origins of non-locality in 
these two approaches. While, non-locality in the many- 
body approach arises from the use of vdW potential with 
a repulsive core and the hypercentrifugal repulsion, that 
in the non local GP approach is built in the choice of the 
two-body potential. 

Although the presence of an intermediate state has al- 
ready been established in earlier works where GP mean- 
field theory has been employed with non-local interac- 
tion [IJl, the effect of the short-range behavior of real- 
istic interatomic potential on the intermediate state has 
not been investigated as yet. Especially in the collapse 
regime, where the atomic cloud becomes highly corre- 
lated, the use of a correlated many-body approach is 
required for the accurate determination of the collapse 
point and various properties of the atomic cluster. The 
non-local interaction in the effective interaction can pre- 
vent a pathological collapse, which results in the GP ap- 
proach with a purely local attractive interaction. The 
comparison of these two approaches is elaborately dis- 
cussed in the present work. In the usual trap setup we 
observe a discontinuous jump between the low-density 
and high-density phases. However, the use of tight traps 
is also interesting and we observe how the discontinu- 
ous jump can be avoided by using microtraps. For large 
number of atoms, the system remains in a quantum self- 
trapped state. 
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